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Abstract 

This  paper  describes  the  application  of  high 
performance  computing  to  the  prediction  of  the  rate 
constants  of  reactions  occurring  in  the  troposphere 
involving  toxic  industrial  compounds.  The  methods  we 
employ  use  a combination  of  quantum  chemistry  and 
quantum  dynamics  to  calculate  the  kinetics  of  the 
reactions  under  investigation.  Our  accomplishments  from 
the  past  year  are  presented  and  discussed. 

1.  Introduction 

The  ability  to  accurately  and  confidently  compute  the 
kinetics  of  potentially  significant  atmospheric  reactions  of 
toxic  industrial  compounds  (TICs)  is  a critical  component 
of  Department  of  Defense  (DoD)  chem/bio  defense 
programs  as  well  as  Homeland  security.  However,  much 
of  the  critical  kinetic  data  (e.g.,  rate  constants  for  target 
reactions)  does  not  exist  and  is  difficult  and  expensive  to 
obtain  experimentally.  Our  unique  approach  is  to  use 
state-of-the-art  computational  quantum  chemistry/ 
quantum  dynamics  models  to  compute  rate  constants 
which  can  then  be  directly  input  into  atmospheric 
chemistry  modules  within  DoD  transport  and  dispersion 
(T&D)  models.  To  this  end,  T&D  models,  such  as 
SCIPUFF  and  the  recently  parallelized  ChemCODE/ 
ChemCONC  (Common  High  Performance  Computing 
Software  Support  Initiative  project  CBD-03)  have  been 
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updated  with  the  ability  to  model  atmospheric  degradation 
of  released  toxic  chemicals.  This  addition  has  resulted  in 
the  improved  ability  of  the  models  to  handle  the  fate  and 
interaction  of  TICs  in  the  atmosphere,  resulting  in  a more 
realistic  representation  of  the  contamination  zone.  The 
methods  we  apply  can  be  used  not  only  to  predict  the 
changing  health  threat  in  both  an  urban  and  battlefield 
environment,  but  can  also  predict  the  nature  of  the 
resulting  byproducts.  This  information  would  be  valuable 
in  programming  advanced  sensors  to  warn  of  the  release 
of  a toxic  compound  even  if  the  material  has  degraded  by 
the  time  it  reached  the  detection  array.  In  addition,  the 
results  of  these  calculations  can  also  be  used  to  make  a 
more  scientifically  defensible  selection  of  simulants  for 
challenging  detectors,  as  well  as  assisting  in  the 
evaluation  of  both  individual  and  collective  protection 
systems.  These  calculations  will  also  be  used  to  provide  a 
sound  theoretical  underpinning  to  the  development  of 
physical  properties  data,  which  is  critical  for  choosing  the 
proper  simulant  for  use  in  synthesis,  fate,  and 
decontamination  studies,  as  well  as  in  the  assessment  of 
new  potential  threat  agents  and  the  selection  of  improved 
decontamination  concepts. 

2.  Methodology 

A team  of  approximately  seven  PhD  level 
computational  scientists  at  the  US  Air  Force,  US  Army 
Research  Laboratory  (ARL),  ENSCO,  Inc.  in  Melbourne, 
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FL,  and  BWD  Associates  in  Gainesville,  FL,  has 
demonstrated  that  accurate  rate  constants  can  be 
calculated  using  a combination  quantum 

chemistry/quantum  dynamics  (QC/QD)  approach.  This 
method  is  computationally  intensive;  it  requires  the  use  of 
high  level  ab  initio  quantum  chemistry  to  calculate 
accurate  structures  of  molecules  at  the  various  important 
stationary  points  along  the  minimum  energy  path  that 
connects  reactants  to  products.  The  computed  energies  of 
the  various  stationary  points  are  also  critically  important. 
An  error  of  ±0.4  kcal/mol  in  the  computed  kinetic  barrier 
height  (the  difference  in  the  computed  energies  of  the 
saddle  point  and  reactant)  can  lead  to  a doubling  or 
halving  of  the  computed  rate  at  298K.  This  clearly  shows 
that  sub-kcal/mol  accuracy  is  required  of  this  energy 
difference,  which  is  by  current  QC  standards,  a very 
challenging  constraint.  To  achieve  this  level  of  accuracy, 
we  use  an  energy  extrapolation  which  accounts  for  the 
well  understood  and  predictable  deficiencies  in  the  QC 
methods  used. 

Additional  information  characterizing  the  transition 
state  is  also  needed;  i.e.,  the  negative  (repulsive)  CCSD 
curvature  of  the  minimum  energy  path  at  the  transition 
state.  Since  a lull  CCSD  vibrational  frequency 
calculation  can  not  be  completed  in  a timely  manner,  we 
compute  the  CCSD  curvature  of  the  repulsive  mode  at  the 
transition  state  by  fitting  a harmonic  potential  to  three 
CCSD  single  point  energies  (one  at  the  transition  state  and 
one  on  each  “side”  of  the  transition  state).  In  addition  to 
the  curvature  of  the  repulsive  mode,  a description  of  the 
soft  modes,  in  the  neighborhood  of  the  transition  state, 
must  also  be  known.  Once  again,  the  molecular  size  (i.e., 
number  of  atoms  and  number  of  basis  functions)  of  the 
compounds  under  study  precludes  the  use  of  coupled 
cluster  theory  to  calculate  this  information.  Instead,  we 
use  a neglect  of  diatomic  differential  overlap  (NDDO) 
Hamiltonian  that  is  reparameterized  to  accurately 
reproduce  coupled  cluster  theory.  Specifically,  given  the 
transition  state  geometry  computed  at  the  coupled  cluster 
level  and  the  curvature  of  the  repulsive  mode  computed 
using  the  three  point  fit  method,  we  vary  the  parameters  of 
the  AMI  variant  of  the  NDDO  Hamiltonian  using  a 
combination  of  simulated  annealing  and  genetic 
algorithms  such  that  the  transition  state  geometry  and  the 
vibrational  modes  generated  using  the  reparameterized 
NDDO  Hamiltonian  compare  with  those  from  coupled 
cluster  theory.  Finally,  this  information  is  provided  as 
input  to  a newly  developed  Semi-Classical  Flux-Flux 
Autocorrelation  Function  (SCFFAF)  approach  to  chemical 
dynamics  which  computes  the  kinetic  rate  curve  of  the 
reaction  under  investigation  over  a wide  temperature 
range.  Details  of  this  methodology  can  be  found  in 
Reference  1.  In  summary,  the  major  steps  required  to 
compute  a rate  constant  are  as  follows: 

1 .  Identify  compound  and  reaction  to  study. 


2.  Optimize  reactant  structures  at  MP2/6-31++G** 
level  of  theory  to  determine  which  conformations 
are  energetically  accessible  at  ambient 
temperature. 

3.  Using  chemical  intuition  and  the  parent 
structures  identified  above,  we  generate  guess 
structures  of  likely  transition  states  and  optimize 
at  MP2/6-3 1 1++G**.  This  typically  can  result  in 
up  to  as  many  as  30  potential  transition  states  of 
interest. 

4.  Once  the  transition  states  structures  have 

converged,  the  activation  energy  (energy  of  the 
transition  state  minus  energy  of  reactants)  for 
each  converged  transition  state  is  determined. 
For  those  transition  states  that  have  an  activation 
energy  of  <3  kcal/mol,  we  initiate  a MP2/6- 
31++G**  vibrational  frequency  calculation.  This 
calculation  serves  a two-fold  purpose:  a)  it 

verifies  that  the  computed  structure  is  indeed  a 
first-order  transition  state  (i.e.,  vibrational 
analysis  contains  one  and  only  one  repulsive 
mode)  and  b)  generates  a guess  Hessian  for  the 
next  series  of  calculations. 

5.  Optimize  selected  transition  states  and  reactant 
structures  at  CCSD/6-31++G**. 

6.  Run  the  energy  extrapolation  calculations  on  the 
converged  transition  states  and  reactants. 

The  extrapolation  method  consists  of  the  following: 

• Single  point  energy:  QCISD(T)/6-3 1 lg** 

(QCI) 

• Single  point  energy:  MP2/6-311g**  (SML) 

• Single  point  energy:  MP2/6-311+G(3df,2p) 

(LRG) 

• Zero  Point  Energy  (ZPE)  from  a vibrational 
frequency  analysis:  MP2/6-31++G** 

The  extrapolated  energy  = E(QCI)  + E(LRG)  - 
E(SML)  + ZPE 

7.  Compute  the  CCSD  curvature  of  the  transition 
state  using  the  three  point  fit  methodology. 

8.  Use  the  CCSD  structure  and  curvature 
information  to  generate  the  Transfer 
Hamiltonian. 

9.  Run  the  Transfer  Hamiltonian  calculations  to 
characterize  the  region  around  the  CCSD/6- 
31++G**  transition  state. 

10.  Compute  the  kinetic  rate  curve  using  the 
SCFFAF  methodology. 

3.  Results 

During  FY  07,  we  have  calculated  temperature 
dependent  rates  for  the  following  hydrogen  abstraction 
reactions: 
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• Reaction  1:  Dimethyl  phosphonate  [DMHP, 

(CH30)2P(0)H]  + OH  radical 

• Reaction  2:  Dimethyl  methylphosphonate 

[DMMP,  (CH30)2P(0)CH3]  + OH  radical 

• Reaction  3:  Diethyl  methylphosphonate  [DEMP, 
(CH3CH20)2P(0)CH3]  + OH  radical 

In  this  manuscript  we  summarize  the  results  of  the 
DMHP  reaction121. 

3.1.  Reaction  1:  DMHP  + OH 

When  released  into  the  troposphere,  dimethyl 
phosponate  (DMHP)  will  degrade,  for  example,  by 
reacting  with  hydroxyl  (OH)  radicals.  The  products  of 
this  hydrogen  abstraction  reaction  are  a radical  DMHP 
species  and  water.  Our  calculations  found  four  transition 
states  that  are  accessible  at  ambient  temperatures  and  are 
thus  the  major  contributors  to  the  total  rate.  Figures  1—3 
show  three  different  transition  states  for  the  abstraction  of 
a hydrogen  atom  from  one  of  the  methyl  groups  on 
DMHP.  Figure  4 shows  the  transition  state  for  the 
abstraction  of  the  hydrogen  atom  that  is  bonded  to  the 
phosphorus  atom.  Figure  5 shows  the  computed 
temperature  dependent  rate  constants  for  each  of  these 
mechanisms,  as  well  as  the  total  rate,  which  is  simply  the 
sum  of  the  individual  rates.  Our  calculated  total  rate 
agrees  quite  favorably  with  the  experimental  rate131, 
shown  in  Figure  5.  Note  that  each  of  the  experimental 
data  points  is  an  individual  measurement  at  a single 
temperature,  whereas  our  calculations  give  the  full 
temperature  dependent  rate. 

4.  Conclusions 

When  released  into  the  troposphere,  volatile  organic 
compounds,  such  as  the  organophosphorus  compounds 
studied  here,  typically  react  with  OH  radicals.  In  this 
work,  we  have  used  a combination  of  quantum  chemistry 
and  quantum  dynamics  to  theoretically  predict  the 
temperature  dependent  rates  for  the  abstraction  of 
hydrogen  atom  atoms  by  OH  radicals  for  a series  of 
organophosphonates.  These  computed  rate  constants  will 
be  used  as  input  into  atmospheric  chemistry  modules 
within  DoD  T&D  models,  such  as  SCIPUFF  and 
ChemCODE/ChemCONC.  The  addition  of  finite 
chemistry  will  result  in  the  improved  ability  of  the  models 
to  handle  the  fate  and  interaction  of  TICs  in  the 
atmosphere,  resulting  in  a more  realistic  representation  of 
the  contamination  zone. 
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Figure  1.  DMHP  + OH  TS1 


Figure  3.  DMHP  + OH  TS2p 
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Figure  4.  DMHP  + OH  TS3 


Temperature  pq 

Figure  5.  Temperature  dependent  rate  constants  for 
the  reaction  of  DMHP  + OH 
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